Demographic DC tidycensus data analysis

Install and load packages

Load required packages for the specific analysis.

library(tidycensus) # Access census data
library(tidyverse)
library(sf)         # Working with simple features - geospatial features
library(tmap)
library(lubridate)  
library(janitor)    # To clean names of variables in the data
library(here)       # A package I haven't taught you about before that doesn't do much, but ....

Going to need an API key to access the data. We can acquire a key from the following website: http://api.census.gov/data/key_signup.html

census_api_key("ef8a7b63162e3c4881110f495943872be80639fd", overwrite=TRUE)

To install your API key for use in future sessions, run this function with `install = TRUE`.

Required Census codes for demographic variables. Selection is based on the most popular demographic values: Age, Gender, Ethnicity, Income

• B01002_001 DC Age Distribution
• B01001_002 Sex (M)
• B01001_026 Sex (F)
• B01001_001 Total population
• B01001H_001 Population (white)
• B02009_001 Population (black)
• B19013_001 Median household income distribution in D.C.

cens_vars <- load_variables(2018, "acs5", cache = TRUE)

Getting together the 2018 DC Census Data

df_census <- get_acs(geography = "tract",                  
                  variables=c("median_inc"="B06011_001",
                              "pop_tot"="B01001_001",
                              "pop_white"="B01001H_001",
                              "pop_black"="B02009_001",
                              "age"="B01002_001",
                              "male"="B01001_002",
                              "female"="B01001_026"),
                  state="DC",geometry=TRUE,year=2018) %>%
  clean_names()

Getting data from the 2014-2018 5-year ACS
Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
Using FIPS code '11' for state 'DC'

Now we should classify and plot the data frame we just created

class(df_census)

[1] "sf"         "data.frame"

plot(df_census)

Now we should clean the census data frame and add percentages for population and sex variables.

df_cens=df_census %>% select(-moe) %>% spread(variable,estimate) %>%
  mutate(pop_black_pct = pop_black/pop_tot, 
         pop_white_pct=pop_white/pop_tot, 
         male_pct=male/pop_tot,
         female_pct=female/pop_tot)

Visualize demographic outcome in maps.

# Visualize median income distribution in Washington D.C.
tm_shape(df_cens) +tm_polygons("median_inc",alpha=.5)

# Visualize age distribution in Washington D.C.
tm_shape(df_cens) +tm_polygons("age",alpha=.5)

# Visualize population distribution in Washington D.C.
tm_shape(df_cens) +tm_polygons("pop_tot",alpha=.5)

# Visualize black population distribution in Washington D.C.
tm_shape(df_cens) +tm_polygons("pop_black",alpha=.5)

# Visualize white population distribution in Washington D.C.
tm_shape(df_cens) +tm_polygons("pop_white",alpha=.5)

# Visualize male population distribution in Washington D.C.
tm_shape(df_cens) +tm_polygons("male",alpha=.5)

# Visualize female distribution in Washington D.C.
tm_shape(df_cens) +tm_polygons("female",alpha=.5)

Now we need to read in and classify the DC neighborhood data

neigh=st_read(here("data", "raw", "DC_Health_Planning_Neighborhoods.geojson")) %>% clean_names()

Reading layer `DC_Health_Planning_Neighborhoods' from data source 
  `C:\Users\Tarak\Documents\CollegeStuff\SeniorYear\Fall2021\DS241\BikeProject\data\raw\DC_Health_Planning_Neighborhoods.geojson' 
  using driver `GeoJSON'
Simple feature collection with 51 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.11976 ymin: 38.79165 xmax: -76.9094 ymax: 38.99556
Geodetic CRS:  WGS 84

class(neigh)

[1] "sf"         "data.frame"

Now reading in the COVID 19 DC data

df_c=readxl::read_excel(here("data","raw","neigh_cases.xlsx"),
                            col_types = c("date", "text", "numeric")) %>% 
  clean_names()

# We should filter the data so we can better analyze it
df_cases=df_c %>%
  filter(as_date(date) == "2021-11-17") %>% 
  separate(neighborhood,into=c("code","name"),sep = ":") %>%
  mutate(code=case_when(
    code=="N35" ~"N0",
    TRUE ~ code
  ))

Now we progress to joining these data frames to accompish the project task

# Need to join neighborhood data with covid case data
neigh2=left_join(neigh,df_cases,by=c("code")) 
tmap_mode("view")

tmap mode set to interactive viewing

# Create map of total positive covid cases in DC
tm_shape(neigh2) +tm_polygons("total_positives", alpha=.5)

  tm_shape(neigh2) +tm_borders(col="blue",lwd=5,alpha=.2)+
  tm_shape(df_cens) +tm_borders(col="red",lwd=1,alpha=.3)

#df_j=st_join(df_cens,neigh2,prepared=FALSE)

df_cens_adj=df_cens %>% st_transform(4326)

df_j=st_join(df_cens_adj,neigh2,largest=TRUE)

Warning: attribute variables are assumed to be spatially constant throughout all geometries

df_final=df_j %>% select(median_inc, pop_tot, pop_black, pop_white, age, male, female, total_positives, objectid) %>%
  group_by(objectid) %>%
  summarise(pop_n=sum(pop_tot),
            pop_black_n=sum(pop_black),
            pop_white_n=sum(pop_white),
            male_n=sum(male),
            female_n=sum(female),
            age_n=sum(age)/n(),
            adj_median_income=sum(pop_tot*median_inc)/pop_n,
            tot_positives_n = sum(total_positives)/n())
            
plot(df_final)

tm_shape(df_final)+tm_polygons(c("adj_median_income","tot_positives_n", "pop_black_n", "pop_white_n"))